function [parms_fitted, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_Nstate(parms, kappa_guess, b_guess, eta_guess, sigma_z_guess, chi_guess, x_bar_guess, log_sigma_x_guess, ...
        U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_vol_target, z_phi_pct, E_phi_z_phi_pct, ave_mu3_vol_target, NumTries, fsolve_opts)
% target moments:
% E[U], vol[u], vol(N*w)/vol(Y), vol(dlogY), vol(ylds), vol(mu3_1y)
% parameters:
% kappa, b, eta, z(1), chi, sigma_x_REC
% sigma_x(N,z) = exp(a + sigma_x_REC*1{z==1})

    if sigma_z_guess<=0; error('sigma_z_guess should be >0'); end
    if eta_guess<=0 || eta_guess>=1; error('check eta_guess'); end
    if kappa_guess<=0; error('check kappa0_guess'); end
    if b_guess<=0; error('check b_guess'); end
    if chi_guess>1 || chi_guess<0; error('check chi_guess'); end

%     FLAG_DEBUG_LOG = true;
    FLAG_DEBUG_LOG = false;
    
%     fixedpt_tol = 1e-8;
    % fixedpt_tol = 1e-7;
    fixedpt_tol = 1e-6; % speed things up

    inflation_process.mu_pi = 6.1592e-04;
    inflation_process.rho_pi = 0.81;
    inflation_process.sigma_pi = 0.00245;
    inflation_process.theta_pi = -0.338;
    inflation_process.loading_dlogC = -0.035;

    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;
    [NN_mesh,xx_mesh] = meshgrid(parms.grid_N(:),parms.grid_x(:));
    Pzz_Nxzz = permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]);
    mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
    var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi + inflation_process.theta_pi^2)...
        *inflation_process.sigma_pi^2/(1 - inflation_process.rho_pi^2);
    var_e = inflation_process.sigma_pi^2;
    cov_Xe = inflation_process.sigma_pi^2;
    
    N_Nxz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]);
    N_Nz = repmat(parms.grid_N(:),[1 parms.Nz]);
    U_Nz = 1 - N_Nz;
    U_Nxz = 1 - N_Nxz;
    x_Nxz = permute(repmat(parms.grid_x(:),[1 parms.NN parms.Nz]),[2 1 3]);
    
    phi_of_x = @(x) 1./(1 + exp(-x));
    phi_Nxz = phi_of_x(x_Nxz);
    log_phi_Nxz = log(phi_Nxz);
    input_mu3_Nxz = (1 - U_Nxz).*U_Nxz.*(1 - 2*U_Nxz).*(log_phi_Nxz.^3);
    
    FLAG_SUCCESSFUL_FSOLVE = false;
    
    for iter = 1:NumTries
        fprintf('iter %g with kappa_guess=%g, b_guess=%g, eta_guess=%g, sigma_z_guess=%g, chi_guess=%g, sigma_x_guess=%g, x_bar_guess=%g.\n', ...
            iter, kappa_guess, b_guess, eta_guess, sigma_z_guess, chi_guess, exp(log_sigma_x_guess), x_bar_guess);
        [XFit, FitVal, EXITFLAG] = fsolve(@FitFun,[log(kappa_guess); ...
            log(b_guess); log(1/eta_guess-1); log(sigma_z_guess); ...
            log(1/chi_guess - 1); log_sigma_x_guess; x_bar_guess],fsolve_opts);
        if EXITFLAG>0
            fprintf('Found roots on iteration %g.\n', iter);
            FLAG_SUCCESSFUL_FSOLVE = true;
            break;
        else
            fprintf('iter %g, fitted vals (not solution):\n kappa=%g, b=%g, eta=%g, sigma_z=%g, chi=%g, sigma_x=%g, x_bar=%g.\n', ...
                iter, exp(XFit(1)), exp(XFit(2)), 1/(1 + exp(XFit(3))), ...
                exp(XFit(4)), 1/(1 + exp(XFit(5))), exp(XFit(6)), XFit(7));
            kappa_guess = 0.1 + 5*rand;
            b_guess = 0.3 + 0.6*rand;
            eta_guess = 0.1*rand;
            sigma_z_guess = rand*0.05;
            chi_guess = 1/(1+4*rand);
            log_sigma_x_guess = log(0.01 + 0.2*rand);
            x_bar_guess = -log(1/(0.85 + 0.1*rand) - 1);
        end
    end
    
    % do a last compute
    parms_fitted = parms;
    parms_fitted.kappa(:) = exp(XFit(1));
    parms_fitted.b(:) = exp(XFit(2));
    parms_fitted.eta(:) = 1/(1 + exp(XFit(3)));
    parms_fitted.sigma_z = exp(XFit(4));
    parms_fitted.chi = 1/(1 + exp(XFit(5)));
    parms_fitted.sigma_x(:) = exp(XFit(6));
    parms_fitted.x_bar = XFit(7);

    [parms_fitted.grid_z,parms_fitted.StateTransitionProbs] = ...
        rouwenhorst(0,parms_fitted.rho_z,parms_fitted.sigma_z^2,parms_fitted.Nz);

    fprintf('kappa=%g, b=%g, eta=%g, sigma_z=%g, chi=%g, sigma_x=%g, x_bar.\n', ...
        parms_fitted.kappa(1), parms_fitted.b(1), parms_fitted.eta(1), ...
        parms_fitted.sigma_z, parms_fitted.chi, exp(XFit(6)), XFit(7));

    gesol = SolveFixedPoint_Model_20220926(parms_fitted,fixedpt_tol);
    parms_fitted.FitVal = FitVal;

    if ~FLAG_SUCCESSFUL_FSOLVE
        warning('Model fit: not all moments are satisfied.');
    end

    function Z = FitFun(X)
        
        Z = NaN*zeros(7,1);
        
        parms_temp = parms;
        parms_temp.kappa(:) = exp(X(1));
        parms_temp.b(:) = exp(X(2));
        parms_temp.eta(:) = 1/(1 + exp(X(3)));
        parms_temp.sigma_z = exp(X(4));
        parms_temp.chi = 1/(1 + exp(X(5))); 
        parms_temp.sigma_x(:) = exp(X(6));
        parms_temp.x_bar = X(7);

        [parms_temp.grid_z,parms_temp.StateTransitionProbs] = rouwenhorst(0,parms_temp.rho_z,parms_temp.sigma_z^2,parms_temp.Nz);
        
        z_Nxz = permute(repmat(parms_temp.grid_z(:),[1 parms.NN parms.Nx]),[2 3 1]); % note: this must be updated when z changes!!
        Y_Nxz = exp(z_Nxz).*N_Nxz;

        gesol_temp = SolveFixedPoint_Model_20220926(parms_temp,fixedpt_tol);
        if gesol_temp.VFIerr>fixedpt_tol
            warning('FitModel20220926_ver1: VFI error = %g is greater than tolerance of %g.', ...
                gesol_temp.VFIerr, fixedpt_tol)
            bComputeMoments = false;
        else
            bComputeMoments = true;
        end
        
        %% Newer version
        
        if bComputeMoments

        [Phi_Nxz,ss_dist_Nxz] = MakeSparseTransMatrix_Nxz(parms_temp,gesol_temp);

        % NEW CODE: quarterly average
        [mean_U, std_U] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,U_Nxz);
        [mean_Y, std_Y] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,Y_Nxz);
        [mean_Nw, std_Nw] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,N_Nxz.*gesol_temp.wages);
        ratio_stdNw_stdY = std_Nw/std_Y;
        
        % quarterly output vol (exact), from library\exact moments
        vol_gy = moment_quarterly_output_growth_vol_Nxz(parms_temp,gesol_temp,Phi_Nxz,ss_dist_Nxz);
        
        % cross-sectional mu3 vol
        std_mu3_cons_1y = CalcTimeSeriesDiffVol(input_mu3_Nxz,ss_dist_Nxz,Phi_Nxz,12);

        % nominal rate
        x_next_Nxzz = gen_law_of_motion_x(parms_temp);
        Nnext_Nxz = gen_law_of_motion_N(parms_temp,gesol_temp);

        C_next = zeros(NN,Nx,Nz,Nz);
        for iznext = 1:parms.Nz
            C_next(:,:,:,iznext) = reshape(interp2(NN_mesh,xx_mesh,gesol_temp.C(:,:,iznext)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iznext)),[NN Nx Nz]);
        end
        inflation_process_temp = inflation_process;
        inflation_process_temp.F_Nxzz = inflation_process_temp.loading_dlogC*(C_next - repmat(sum(Pzz_Nxzz.*C_next,4),[1 1 1 Nz]));
        
        gesol_temp.nom_bonds = PriceNomBonds_Model20220926(parms_temp, gesol_temp, inflation_process_temp, 12);

        NY = 5;
        [Phi_NxYz,grid_Y,prob_NxYz] = MakeSparseTransMatrix_NxYz(parms_temp,gesol_temp,inflation_process_temp.rho_pi,...
            inflation_process_temp.F_Nxzz,NY,ss_dist_Nxz);

        idx = find(gesol_temp.nom_bonds.maturities==12); mat_yr = 1;

        temp_input = permute(repmat(-log(gesol_temp.nom_bonds.P_A(:,:,:,idx))/mat_yr,[1 1 1 NY]),[1 2 4 3]) ...
            + permute(repmat(gesol_temp.nom_bonds.P_C_a(idx)*grid_Y(:)/mat_yr,[1 NN Nx Nz]),[2 3 1 4]);
        [temp_mean,temp_std] = Calc_monthly_vol(prob_NxYz,temp_input);
        var_Xe_component = (gesol_temp.nom_bonds.P_B_b(idx)^2*var_X + gesol_temp.nom_bonds.P_B_c(idx)^2*var_e ...
            + 2*gesol_temp.nom_bonds.P_B_b(idx)*gesol_temp.nom_bonds.P_B_c(idx)*cov_Xe)/(mat_yr^2);

        nom_yld_std = 100*sqrt(temp_std^2 + var_Xe_component);


        prob_x_z = squeeze(sum(ss_dist_Nxz,1));
%         prob_z = sum(prob_x_z,1);
% 
%         % compute E[phi(t)|z(t)>=z] and Prob(z(t)>=z)
%         E_phi_given_geq_z = zeros(parms_temp.Nz,1);
%         P_gez_z = zeros(parms_temp.Nz,1);
%         for iz = 1:parms_temp.Nz
%             P_gez_z(iz) = sum(prob_z(iz:end));
%             E_phi_given_geq_z(iz) = sum(prob_x_z(:,iz:end).*repmat(phi_of_x(parms_temp.grid_x(:)),[1 (parms_temp.Nz-iz+1)]),'all')...
%                 /P_gez_z(iz);
%         end
% 
%         if 1-parms.RECESSION_PROB_THRESH<min(P_gez_z) || 1-parms.RECESSION_PROB_THRESH>max(P_gez_z)
%             error('Increase Nz.');
%         end
%         E_phi_boom = interp1(P_gez_z,E_phi_given_geq_z,1-parms.RECESSION_PROB_THRESH);

        z_phi = norminv(z_phi_pct/100,0,parms_temp.sigma_z/sqrt(1-parms_temp.rho_z^2)); % AR(1) process
        E_phi_given_z = zeros(parms_temp.Nz,1);
        
        for iz = 1:parms_temp.Nz
            E_phi_given_z(iz) = sum(prob_x_z(:,iz).*phi_of_x(parms_temp.grid_x(:)))/sum(prob_x_z(:,iz));
        end
        
        E_phi_boom = interp1(parms_temp.grid_z,E_phi_given_z,z_phi);
        
        Z(1) = (mean_U - U_mean_target)/U_mean_target;
        Z(2) = (std_U - U_std_target)/U_std_target;
        Z(3) = (ratio_stdNw_stdY - wagebillratio_target)/wagebillratio_target;
        Z(4) = (vol_gy - vol_gy_target)/vol_gy_target;
        Z(5) = (nom_yld_std - nom_yld_1y_vol_target)/nom_yld_1y_vol_target;
        Z(6) = (std_mu3_cons_1y - ave_mu3_vol_target)/ave_mu3_vol_target;
        Z(7) = (E_phi_boom - E_phi_z_phi_pct)/E_phi_z_phi_pct;

        if FLAG_DEBUG_LOG
            fprintf('\nkappa=%4.10f,\nb=%4.10f,\neta=%4.10f,\nsigma_z=%4.10f,\nchi=%4.10f,\nsigma_x=%4.10f,\nx_bar=%4.10f.\n',...
                parms_temp.kappa(1),parms_temp.b(1),parms_temp.eta,parms_temp.sigma_z,parms_temp.chi(1),parms_temp.sigma_x(1),...
                parms_temp.x_bar(1));
            fprintf('\nE[U]=%g,\nstd[U]=%g,\nw_elas=%g,\nstd[gY]=%g,\nstd[yld]=%g,\nstd[mu3]=%g,\nE[phi|boom]=%g.\n\n',...
                mean_U, std_U, ratio_stdNw_stdY, vol_gy, nom_yld_std, std_mu3_cons_1y, E_phi_boom);
        end
        
        end
        
    end

end